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We review the theory of optimal polynomial and rational Chebyshev approximations, and Zolotarev's formula 
for the sign function over the range e < \z\ < 1. We explain how rational approximations can be applied to 
large sparse matrices efficiently by making use of partial fraction expansions and multi-shift Krylov space solvers. 



1. Introduction 

There are many situations in which it is desir- 
able to evaluate a function of a matrix. For in- 
stance, in lattice quantum held theory it is some- 
times desirable to evaluate the square root of a 
discretised Dirac operator Fj in order to calcu- 
late the effects of varying the number of fcrmionic 
flavours [1I2I3I4I5] , or to construct a good approx- 
imation to Neuberger's operator for Ginsparg- 
Wilson fermions by evaluating the sign function 
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2. Matrix Functions 

As a general mathematical problem we need to 
dehne what we mean by the generalisation of a 
scalar function to a corresponding matrix-valued 
function on matrices. This can be done in a va- 
riety of ways |10| . but for our purposes it suf- 
fices to restrict our attention to the special case 
of an Hermitian matrix H, which can be trans- 
formed into real diagonal form by a unitary trans- 
formation 1 H = UDU^ . Defining a function of a 
real diagonal matrix to be a diagonal matrix in 
the obvious way, f{D)a = f(Da), we may define 
f(H) ee Uf(D)W. 

We next address the question of how to com- 
pute good numerical approximations to f(H) 
cheaply for a useful class of functions. One way 
of doing this is to find a polynomial or rational 
function which is a good approximation over a 
compact interval R containing the spectrum of 

1 Although the lattice Dirac operator is not hermitian 
the assumption of 75 hermiticity implies that 075 is, and 
this may be used instead. 



H . This seems a promising approach because it 
is well known that rational approximations are an 
effective way of computing scalar functions on a 
digital computer 111112113114) . 

The reason why matrix polynomials are cheap 
to evaluate is that they do not require explicit 
diagonalisation of the matrix. From the identities 

U(aD + f3D')U^ = alJDU^ + f3UD'U\ 
UDD'U^ = UDU^UD'U\ 

we see that p(H) = p(UDU r ) = Up(D)W for 
any polynomial p. Similarly, if we observe that 
UD^W = [UDW]- 1 then we see that r{H) = 
r(UDW) = Ur(D)W for any rational function 
r, which is cheap to evaluate if we count matrix 
inversion as a "cheap" operation. 

It is reasonable to ask why we should consider 
matrix inversion as "cheap." There are three 
practical reasons for this 

1. In applications we do not need a function of 
a matrix or the inverse of a matrix per se, 
but only the effect of applying it to some 
vector. Therefore we only need to solve 
a system of linear equations, which is far 
cheaper than finding the full inverse. 

2. There are very efficient methods for solv- 
ing large systems of sparse linear equations. 
We may take advantage of all the work that 
has been done in developing and optimising 
such Krylov space algorithms. 

3. We can expand a rational function in par- 
tial fractions |7I8I15| . and then take advan- 
tage of the structure of Krylov space solvers 
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to compute all the terms in the partial frac- 
tion expansion using the same Krylov space 
|lbll7| . For this to work it is necessary that 
the partial fraction expansion is numerically 
stable, i.e., the terms do not suffer large 
cancellations; and that the denominator of 
the rational function is square-free. For rea- 
sons that are not fully understood both of 
these conditions seem to be satisfied in prac- 
tice. 

Observe that the definition of f(H) only re- 
quires that we know the values that / takes on 
the eigenvalues of H , whereas the polynomial or 
rational approximation must be valid over an in- 
terval containing the spectrum. For certain func- 
tions, such as 1/x, we can find cheaper methods: 
for instance, Krylov space methods are usually 
much better for finding an inverse than using a 
polynomial approximation to 1/x |18| . 

3. Approximation by Polynomials 

3.1. Function Norms 

In order to define more precisely what we mean 
by a "good" approximation p to a continuous 
function / over the unit interval [0, 1] it is use- 
ful to introduce a family of norms, 

\\f\\n=(J dxw(x)\f(x)\ n ^j 

where w is a positive weight function and n > 1. 
We may verify that this does indeed define a norm 
on the space of continuous functions over the unit 
interval by the use of the Cauchy-Schwarz in- 
equality. Some cases of especial interest are when 
w(x) = 1 (absolute norm), w(x) = l/|/(x)| (rela- 
tive norm), n = 2 (Euclidean norm), and n — oo 
(Chebyshev norm) for which 

ll/lloo = max w(x)\f(x)\. 

0<x<l 

In this last case if we choose the function p 
from some special class of continuous functions, 
like polynomials or rational functions, then the p 
which minimises the Chebyshev norm of p — / is 
called the minimax approximation, 

\\p - /||oo = min max w(x)\p(x) - /(x)|. 

p 0<:r<l 1 1 



If we are using floating-point numbers to rep- 
resent components of spinor fields on a com- 
puter then it is most consistent to use a min- 
imax approximation with a relative weight for 
matrix function evaluation. In particular, using 
an approximation with a minimax error of one 
ulp (unit of least precision) means that the errors 
caused by function approximation have no reason 
to be worse than those we already accept by using 
floating-point arithmetic. 

3.2. Weierstrass' theorem 

The fundamental theorem on the approxima- 
tion of continuous functions by polynomials is 
due to Weierstrass, who proved that any continu- 
ous function can be arbitrarily well approximated 
over the unit interval by a polynomial. This the- 
orem is of significance in functional analysis, as it 
shows that the polynomials are dense in the space 
of continuous functions with respect to the topol- 
ogy induced by the Chebyshev norm. An elegant 
proof of Weierstrass' theorem may be obtained by 
the use of Bernstein Polynomials, 

k=Q V ' X 7 

It is simple to show that linin^oo \p n — f\\oo — 0. 

4. Optimal Approximation 

Weierstrass' theorem tells us that for any spec- 
ified tolerance e > we can find a degree n 
such that there is a polynomial p n satisfying 
\\Pn — f\\oo < £• For our purposes the con- 
verse question is of greater importance, namely 
for a given degree n what is the smallest error 
\\Pn — /II oo that can be achieved. This ques- 
tion was addressed by Chebyshev, who provided 
a suprisingly simple answer. For pedagogical rea- 
sons we shall first discuss Chebyshev's theorem 
for the case of approximation by polynomials, 
and then generalise the result to rational func- 
tions (ratios of polynomials); although logically 
the polynomial case is just a special case of the 
rational one where the degree of the denominator 
polynomial happens to be zero. 
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Figure 1. The solid line shows the error e for a 
polynomial approximation p of degree d to the 
continuous function /, which has less than d + 2 
extrema of equal magnitude. This defines a "gap" 
A between the extrema. We may construct the 
polynomial q of degree d shown by the dotted line 
that has opposite sign to the error at each of the 
extrema, and whose magnitude A' < A. The 
polynomial p + q is then a better approximation 
than p, as is obvious from its error w(x)[p(x) + 
q{x) — f(x)\ shown by the dashed line. 



4.1. Polynomials 

Chebyshev proved that for any degree d there 
is always a unique polynomial 2 p that minimises 
|| e|| oo = maxo<a;<i |e(aff) [ , where the error is 
e(x) = w(x)\p(x) — f(x)]. Furthermore this poly- 
nomial is characterised by the criterion that the 
error e(x) takes its maximum absolute value at at 
least d + 2 points on the unit interval (which may 
include the end points of the interval), and the 
sign of the error alternates between its successive 
extrema. 

The proof of Chebyshev's theorem is straight- 
forward, and is illustrated in Figures ^ and 
First we prove that Chebyshev's criterion is neces- 
sary: if the error attains its extreme value at fewer 

2 The polynomial may be of degree less than d; for exam- 
ple the best approximation of any degree to a constant 
is just that constant. A less trivial example is that the 
zero polynomial is the best approximation to sin Nirx over 
[0,1], N 6 N, for all d < N. 




Figure 2. The solid curve is the error e p for a 
polynomial p of degree d satisfying Chebyshev's 
criterion. The dotted line shows the error e p i 
for a hypothetical better approximation p' . The 
dashed line is the polynomial p' — p which must 
have at least d + 1 zeros and therefore vanishes 
identically. 



than d + 2 "alternating" points then the approx- 
imation can be improved. Consider a polynomial 
p for which the error e{x) = w(x)[p(x) — f{x)] 
takes its extreme values at fewer than d+2 points 
(shown by the solid curve in Figure ^1; the next 
largest extremum of the error (which may be 
a local extremum or may occur on the bound- 
ary) must be smaller by some non-zero "gap" 
A. As the d + 1 extrema alternate in sign the 
error must pass through zero between each suc- 
cessive extremum, and there are at most d such 
zeros Zi (if there are several zeros between ad- 
jacent extrema we just chose one of them). We 
can construct a polynomial u of degree d (shown 
by the dotted line in the figure) which is zero at 
each of the Zi. This may be written explicitly as 
u(x) = AY[i(x — z^. The constant A may be 
chosen such that the sign of u is opposite to that 
of the error at each extremuum, and its magni- 
tude A' = ||it||oo = maxo< x <i w(x)\u(x)\ is less 
than the gap A. It follows that the polynomial 
p + u is a better approximation to /, as its error 
p(x) + u(x) — f(x) (the dashed line in the figure) 
is everywhere less than e(x). 

Next we prove the sufficiency of Chebyshev's 



4 



criterion: if the error attains its extreme val- 
ues at exactly d + 2 alternating points then 
it is indeed the optimal approximation. To 
show this assume that there was a polynomial 
p' that furnished a better approximation. The 
solid line in Figure shows the error e p {x) = 
w{x)[p{x) — f(x)] for a polynomial p satisfying 
Chebyshev's criterion. The dotted line shows the 
error e p >(x) = w{x)[p' (x) — f(x)] for the hypo- 
thetical better approximation p' . Since p' is a 
better approximation ||e P '||oo < ||e p ||oo, or equiv- 
alently |e P '(xj)| < | e p (a;^ ) | at each of the d + 2 
extrema Xi of e p (x). Appealing to continuity we 
deduce that e p '(zi) — e p (zi) at at least d + 1 
points Zi between the extrema. From this it fol- 
lows that e p '(zi) - e p (zi) = w{zi)\p'{zi) - /(«<)] - 
w(zi)\p(zi) - f(Zi)\ =w(z i )\p'(z i )-p(z i )] = 0, so 
the polynomial p'(x)—p(x) (shown by the dashed 
line) must have a least d + 1 zeros, but as it is of 
degree d this means it is identically zero by the 
fundamental theorem of arithmetic. Hence p' = p 
and the polynomial p is both optimal and unique. 

4.2. Chebyshev Polynomials 

The Chebyshev polynomials 3 T n (cos9) = 
cos n9 provide a simple example of optimal min- 
imax approximations. T n (x) obviously has n + 1 
extrema which alternate in value between — 1 
and 1 for —1 < x < 1; therefore p n (x) = 
x" — 2 1 ~ n T n (x) is the best polynomial approx- 
imation of degree n — 1 with uniform weight to 
the function x n over the interval [—1,1], for the 
error e n {x) = p n (x) — x n = 2 1 ~ n T n (x) satisfies 
Chebyshev's criterion. The magnitude of the er- 
ror is just 1 1 e„ | |oo = 2 1_ ™ = 2e~" ln2 , so we see 
that in this case the error decreases exponentially 
with respect to n. 

If we consider a function / that has a conver- 



3 It is easy to show that these are polynomials, indeed 

L"/2J k 
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The leading coefficient of T n [x), that is the coefficient of 
x", is thus seen to be Y^k^ \ 2k ) = 2 " _1 - 



gent Taylor series expansion over [—1,1], 
^ /«(0) j 

■fW = .L — — x - 

3=0 J ' 

and we reexpress it as an expansion in Chebyshev 
polynomials 4 

oo 

j=o 

then truncating this expansion at order d gives 
an approximation close to the optimal one of this 
degree. It is in fact the optimal polynomial ap- 
proximation to the truncated Taylor expansion 
of degree d + 1, which is not necessarily exactly 
the optimal one. Likewise, the convergence of the 
truncated Chebyshev expansion may converge ex- 
ponentially in d, but it will not if the original 
Taylor series expansion did not converge expo- 
nentially. 

4.3. Rational Functions 

It is perhaps surprising that Chebyshev's ar- 
gument can be easily extended to the case of 
optimal rational approximations as well. The 
statement of Chebyshev's theorem in this case 
is that for any degree (n, d) there is always 
a unique rational function r that minimises 
|| e|| oo = maxo< x <i |e(x)|, where the error is 
e(x) = w(x)[r(x) — f{x)\. For rational approx- 
imations Chebyshev's crtierion is that the error 
e(x) takes its maximum absolute value at at least 
n + d + 2 points on the unit interval (which may 
include the end points of the interval), and the 
sign of the error alternates between its successive 
extrema. 

The proof is similar to that for polynomials 
given above, and we shall therefore just briefly 
consider the salient differences. To prove that 
Chebyshev's criterion is necessary we consider a 
degree (n, d) rational function r = p/q for which 

4 This is straightforward because they satisfy the orthogo- 
nality relation 

/l ^ 
_j Vl - X 2 

(apart from the trivial case where n = m = for which 
there is an additional factor of two). 
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the error e = r— f attains its extreme value HeHoo 
at fewer than n + d + 2 alternating points. The 
denominator q cannot have any zeros on the unit 
interval, as clearly q = 1 would then lead to a 
smaller error; therefore without loss of generality 
we can assume q(x) > 5 > throughout the inter- 
val. Just as before we may construct the polyno- 
mial u of degree n + d that shares the zeros of the 
e. Since gcd(p, q) ^ we can construct polynomi- 
als s and t of degrees d and n such that u = sp+tq 
using the extended Euclidean algorithm. We then 
define the rational function r' = ^+^- which sat- 

q+es 

isfics 



r'-/l = 



P - g _ P + P _ j 
q + es q q 

e(qt+ps) 
q(q + es) 



q{q 



Choosing the constant e small enough that q(x) + 
es(x) > 5' > for x £ [0,1] and eWu]^ < A85' 
the required result follows. 

The proof of sufficiency in the rational case is 
even simpler. If r = p/q is a, degree (n, d) ra- 
tional approximation whose error w(r — f) has 
n + d+2 alternating extrema of equal magnitude, 
and r' — p 1 jq' is a hypothetical better rational 
approximation of the same degree, then the error 
for r' must cross that for r at least n + d+1 times, 



i.e., the rational function 



must have more 



than n + d zeros. Now, in order for 



P 

q 



p_ 

q' 



pq -pq 
qq' 



to vanish on the unit interval the numerator must 
vanish (the denominator certainly cannot); and 
the polynomial pq' — p' q is of degree n + d and 
thus can only have n + d+1 zeros if it vanishes 
identically, hence r — r' . 

5. Remez Algorithm 

Chebyshev's theorem assures us that there is 
a unique 5 optimal rational approximation of any 
degree to any continuous function over the unit 
interval, and it even provides a simple criterion 
for the behaviour of the error in this case, which 



5 In some cases the optimal approximation may actually 
be of lower degree. 



is certainly of practical use in verifying that a 
putative minimax solution is correct. In fact, the 
proof of the theorem is quite constructive, but it 
does not provide a particularly efficient means of 
computing the optimal solution in practice. 

In general, the computation of the optimal 
Chebyshev rational approximation is often per- 
formed using the Remez algorithm |19I14| . To 
describe the Remez algorithm we define a ref- 
erence X as a set of points {xi € [0,1], i = 
1, . . . , n + d + 2} that we will make converge to 
the alternating set at which the error takes its 
extreme values. The Remez algorithm alternates 
two steps: 

1. New r. Keep X fixed, and choose the 
best rational approximation r which passes 
through the points (xj, f{xi) + (— 1)'A) for 
Xi G X. 

2. New X. Keep r fixed, and choose a new 
reference which is the best alternating set 
for r near to X. 

In slightly more detail: 

1. Interpolate an (n, d) rational function 
r = p/q such that the error e(xj) = 
w( Xi )[f lxi) - r{xi)] = (-1) 2 A. Setting 

P( x ) = TJj=oPi xi and q( x ) = Ej=o9i^ J 
with qd = 1 this gives the n + d + 2 equa- 
tions EjLo^L/X^) _ (-i^/X^)] - 

Y^j=oPj x l = f° r tnc n + d + 2 variables 
Pj, qj, and A. Noting that these equations 
are linear in the pj and qj, we write this 
as the matrix equation Mv = where v is 
the vector (p , . . . ,p n , q , . . . , q n ). This has 
a non-trivial solution for v iff detM = 0, 
which is a polynomial in A. For the values 
of A which are real roots of this polynomial 
we can find the coefficients pj and qj, and 
check that the resulting rational function r 
is valid. 

2. Choose a partition of [0, 1] into intervals 
Ii such that Xi € Ii, and choose a new 
reference X 1 = {x- € h : (-l) i e(x' i ) = 
max :re / i (-l) l e(x)}. 

The Remez algorithm has several drawbacks: 
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• It can sometimes be hard to find a suitable 
initial reference: the location of the zeros of 
Chebyshev polynomials are sometimes sug- 
gested as an ad hoc guess. 

• The construction of the interpolating ratio- 
nal function r can fail for some exceptional 
references. 

• Finding the extrema of e(x) in an interval 
Ii can be tricky if the function / oscillates 
rapidly within the interval. As /, and there- 
fore e, can be any continuous function an ef- 
ficient yet robust algorithm for locating the 
extrema is hard to find. 

• The computation often has to be carried 
out in multiple-precision arithmetic: it is 
not unusual to need 20-30 or more decimal 
digits precision. 

• The rate of convergence is quite slow. 

Of course, if bounds on the spectrum of a fam- 
ily of matrices is known a priori then the cost 
of computing the optimal rational approximation 
is not too significant, as it only has to be done 
once for many matrix function evaluations. This 
is the case, for example, in computing roots of 
the Dirac operator for staggered fermions, whose 
spectrum is bounded above by a constant and be- 
low by the mass. On the other hand, in the case 
of Wilson fermions or massless overlap Ginsparg- 
Wilson fermions a useful lower bound on the spec- 
trum does not exist, and one either has to use a 
higher-degree approximation covering a conserva- 
tive range or one has to recompute the approxi- 
mation "on the fly," both of which may be pro- 
hibitively expensive options. 

6. Zolotarev's Theorem 

It is fortunate that in the most interesting 
cases, namely square roots, inverse square roots, 
and the sign function sgn(a;) the coefficients of 
the optimal Chebyshev rational approximation 
are known analytically. This result is due to 
Zolotarev, who was a student of Chebyshev. 

Zolotarev [^Oj gave an explicit expression for 
optimal rational approximations to the sign func- 
tion. To be precise, he gave an expression of the 
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Figure 3. Sketch of the function r(x) of equa- 
tions Q and ©. The four values ±1 ± A are 
indicated by the horizontal lines, and r(x) crosses 
these lines at ±e and ±1. 



form r(x) = xR(x 2 ) where R is a rational func- 
tion and 

max \r(x) — sgn(x)| = A. (1) 

e<|x|<l 

Note that this degree (n + 1, n) approximation 
is as near to a "diagonal" rational function as 
possible, since r must be odd. Furthermore, 
no approximation with lower degree numerator 
could be better, for otherwise it would violate the 
uniqueness part of Chebyshev's theorem. 

The approximation is illustrated in Figure [3] 
It is clear that f(x) = (1 — A 2 )/r(x) is an opti- 
mal approximation to sgn(a;) of degree (n, n + 1), 
which is singular rather than zero at x = 0. If 
there was a better approximation than r with 
higher degree numerator then f would not be 
unique, again contrary to Chebyshev's theorem. 

6.1. Zolotarev Coefficients 

It is amusing that whereas the properties of 
Chebyshev polynomials follow from the fact that 
cos n0 is a polynomial in cos 8, Zolotarev's theo- 
rem follows from an analogous result for elliptic 
functions. 

The Jacobi elliptic function sn is implicitly de- 
fined as the inverse of an elliptic integral, namely 
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It is a doubly periodic analytic function in the 
complex plane with periods uj = 4K(fc) and u>' = 
2?K(/c'), where K is the complete elliptic integral 



K(k) 



dt 



v /(l-t 2 )(l-^ 2 ) 



and fc 2 + fc' 2 = 1. It can be shown |21l22j that 
any elliptic function with these periods must be 
expressible as a rational function of sn(z; k) and 
the related Jacobi elliptic functions cn(z; k) = 
y/1 — sn(z; k) 2 and dn(z; k) — — k 2 sn(z; fc) 2 ; 
and more specifically that any even elliptic func- 
tion with the same periods is a rational function 
of sn(z;fc) 2 . In fact, any elliptic function with 
equivalent periods, in the sense that they generate 
the same period lattice, are rationally expressible 
in this manner; as are elliptic functions with peri- 
ods which divide u> and uj' , since the latter are also 
periodic with the original periods. Such rational 
relationships are known as modular transforma- 
tions. 

Consider in particular an elliptic function with 
periods uj — uj — 4K(fc) and uj' = uj' /n = 
2iK(k')/n. Such a function is easily constructed 
from sn(z; k) by scaling the argument z by some 
factor 1/M and choosing a new parameter A. The 
function sn(z/M; A) has periods 4L = 4K(A) and 
2iL' = 2i K(A') where A 2 + A' 2 = 1, so the periods 
with respect to z are ALM and 2iL'M, which thus 
requires LM = K{k) and L'M = K(k')/n. The 
even elliptic function sn(z/M; A)/ sn(z; k) must 
therefore be expressible as a rational function R 
of sn(z; fc) 2 , sn(z/M; A) = sn(z; k)R (sn(z; fc) 2 ). 

Matching the poles and zeros, and fixing the 
overall normalisation from the behaviour near 
z — 0, the coefficients of the numerator and de- 
nominator polynomials for R may be expressed 
in terms of Jacobi elliptic functions, 



sn(z/M ; A) 
sn(z; fc) 



[n/2] 



l(z;fc) 2 



in 



i(2iAT'm/n;fc) 2 



1 1 



sn(z;fc) 2 



•(2) 



\(2iK'(m-^)/n;k) 



The quantity M is determined by evaluating 
this identity at the half period K(fc) where 
sn(K(fc);fc) = sn(K(fc)/Af; A) = sn(X; A) = 1. 
Likewise, the parameter A is found by evaluating 
the identity at z = K(fc) + iK(k)/n. 
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Figure 4. The fundamental region for sn(z; fc) is 
shown by the dotted box, and that for sn(z/M; A) 
whose period is divided by n = 5 along the imag- 
inary axis is shown by the dashed boxes. The line 
from —K — iK' through —K and K to K + iK' 
shows the contour in the complex z plane along 
which both functions are real, and sn(z; fc) takes 
values from —1 through — e and e to 1 while 
sn(z/M; A) oscillates about it with amplitude A. 



The fundamental region for sn(z; fc) and 
sn(z/M;A) for n = 5 are shown in Figure 0] 
which also shows the contour in the complex z 
plane along which the modular transformation 
gives the desired real rational relationship be- 
tween these two elliptic functions. Some further 
insight is given by Figure[Sl in which the real part 
of sn(z/M;A) is shown over the same region of 
the complex plane. The imaginary part vanishes 
along the contour described in Figure which is 
also indicated on this plot. From this plot it is 
clear how the two periods of the elliptic functions 
correspond to the "step" of the sign function and 
the oscillations of r about it. 



6.2. Approximations to *fx and \j\fx 
We may rewrite equation Q in the form 



max 

e<W<l 



max 

e<|x|<l 



xR(x 2 ) 



R{x 2 ) - -- 



and, upon setting £ = x 2 , we find the relation 



A = max 



• This shows 
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Figure 5. Surface plot showing the real part 
of sn(z/M;A) over the fundamental region of 
sn(z;fc). The imaginary part vanishes along the 
contour indicated. 



that R(x) is an optimal Chebyshev approxima- 
tion to the function 1/y/x with weight \fx over 
the interval [y/e, 1] , or in other words the opti- 
mal Chebyshev approximation with unform reia- 

|fl(«)-i/Vs| 

max v^<«<i 



tive error A 



We can 



rewrite this as A = max v /^ < ^ <1 J —= 

xR(x) is an optimal Chebyshev approximation to 
y/x over [-^/e, 1] with uniform relative error. 



7. Examples 

7.1. A Simple Numerical Example 

It might be useful to look at an actual exam- 
ple of the optimal minimax rational approxima- 
tions we have been discussing. The optimal de- 
gree (3, 3) approximation to 1/v 7 ^ over the inter- 
val [0.003, 1] is 



1 



0.3904603901 



(x+ 2.3475661045) 
(x + 0.4105999719) 



(x + 0. 1058344600) (x + 0.0073063814) 
X (x + 0.0286165446) (x + 0.0012779193)' 

and, as promised, it has a numerically stable par- 
tial fraction expansion 



1 



0.3904603901 



0.1408286237 



0.0511093775 
+ 0.0012779193 + 
0.5964845033 



0.0286165446 



0.4105999719 



Notice that all the terms are positive, and the 
roots of the denominators are all negative. 

7.2. Comparison of Polynomial and Ratio- 
nal Approximations 

We now briefly compare the quality of polyno- 
mial and rational approximations. It is common 
folklore that even though they require division op- 
erations rational approximations are superior to 
polynomial ones, and this seems to be true for ma- 
trix function approximations too, if use is made 
of the partial fraction expansion and multi-shift 
solver tricks. In a sense, as we mentioned before, 
this is because Krylov space solvers are more ef- 
fective than polynomial approximations to 1/x of 
fixed degree as a means of solving systems of lin- 
ear equations. 

A numerical comparison of the minimax er- 
rors for optimal rational and polynomial ap- 
proximating functions to 1/ ^fx over the interval 
[0.00003, 1] (corresponding to a staggered fermion 
mass parameter of to = 0.025) as a function of ap- 
proximation degree is shown in Figure|B] The key 
points to note are: 

• The error of the rational approxima- 
tions falls exponentially with their degree. 
Machine-precision errors of one ulp (about 
10~ 7 for 32-bit IEEE floating-point arith- 
metic and 10 -15 for 64-bit arithmetic) are 
easily achieved for relatively low degree ra- 
tional functions. 

• In order to reach small errors, of the scale of 
the floating-point error scale of one ulp, us- 
ing polynomials requires extemely high de- 
gree polynomials. Considerable care needs 
to be taken with regard to rounding errors 
in these cases. 
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Figure 6. Comparison of minimax errors for opti- 
mal rational and polynomial approximating func- 
tions to x^ 1 !" 1 over the interval [0.00003,1] (corre- 
sponding to a staggered fermion mass parameter 
m = 0.025) as a function of approximation de- 
gree. 



Figure 7. The minimax error A for Zolotarev ap- 
proximations of degree n over the interval [e, 1]. 
The size of the squares indicates the size of the 
rounding errors in evaluating the Zolotarev ratio- 
nal functions in single-precision arithmetic. 



• The amount of work needed to apply the ra- 
tional approximations of a sparse matrix to 
a vector, using a partial fraction expansion 
and multi-shift solver, depends essentially 
on the degree of the rational function times 
the dimension of the Krylov space required. 

The dependence of the minimax error A for the 
Zolotarev rational approximations on the approx- 
imation interval [e, 1] and the degree n is shown in 
Figure [7| This plot also attempts to indicate the 
numerical errors in evaluating the rational func- 
tion in single-precision arithmetic by the size of 
the squares. Again, the rapid convergence and 
stability of the approximations is manifest. The 
errors A appear to be consistent with an asymp- 
totic formula of the form A oc e n ' ln e ^2] ■ 

At the end of section [5] we noted that it was 
expensive to use the Remez algorithm to com- 
pute optimal rational approximations "on the fly" 
in the case where we do not have a useful lower 
bound on the spectrum of a family of Dirac oper- 
ators. For the interesting cases where Zolotarev's 
formula is applicable we can compute the coef- 



ficients in the rational approximation efficiently 
"on the fly" if we measure the smallest eigenvalue 
of each Dirac operator. The method of computing 
the Zolotarev coefficients rapidly and accurately 
making use of Gauss' arithmetico-geometric mean 
is explained in |22| . 



8. Conclusion 

We reviewed the subject of minimax approxi- 
mation theory, and explain why it is particularly 
effective for the approximation of matrix func- 
tions. We demonstrate why rational approxima- 
tions are often superior to polynomial ones, and 
how the combination of partial fraction expan- 
sions and multi-shift Krylov space solvers allows 
them to be applied efficiently to large sparse ma- 
trices. Finally, we have outlined how the optimal 
approximation for the sign function may be found 
in closed form using elliptic functions. 
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